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The Monte-Carlo method is applied to the Polyakov-loop extended Nambu-Jona-Lasinio (PNJL) 
model. This leads beyond the saddle-point approximation in a mean-field calculation and introduces 
fluctuations around the mean fields. We study the impact of fluctuations on the thermodynamics 
of the model, both in the case of pure gauge theory and including two quark flavors. In the two- 
flavor case, we calculate the second-order Taylor expansion coefficients of the thermodynamic grand 
canonical partition function with respect to the quark chemical potential and present a comparison 
with extrapolations from lattice QCD. We show that the introduction of fluctuations produces only 
small changes in the behavior of the order parameters for chiral symmetry restoration and the 
deconfinement transition. On the other hand, we find that fluctuations are necessary in order to 
reproduce lattice data for the flavor non-diagonal quark susceptibilities. Of particular importance 
are pion fields, the contribution of which is strictly zero in the saddle point approximation. 

I. INTRODUCTION 

Understanding the thermodynamics of strongly interacting matter is a persistent challenge. The QCD phase diagram 
in the plane of temperature T and baryonic chemical potential fi features three regions that are accessible with different 
strategies. Along the temperature axis with (i = 0, lattice QCD provides an ab-initio, non-perturbative framework, 
limited only by the available computing power. At asymptotically large chemical potential fj,, perturbative QCD is 
applicable. In the broad range between these extremes, models based on the symmetries and symmetry breaking 
patterns of QCD are useful tools for orientation. 

Different models and approaches have been developed in the last few years for this task [IHB]. A remarkably 
successful model in this context is the Polyakov-loop extended Nambu-Jona-Lasinio (PNJL) model. The NJL model 
[7HS] offers a schematic, but nonetheless quite realistic, picture of the basic dynamics behind spontaneous chiral 
symmetry breaking: the chiral condensate as the order parameter and the pions as Goldstone bosons emerge from 
a chiral invariant, local four-point interaction between quarks. Quark confinement is implemented in addition by 
introducing the Polyakov loop as the order parameter for the confmcmcnt-dcconnnement transition in the pure gauge 
theory and using the minimal gauge invariant coupling to quarks. This produces a dynamical entanglement of the 
chiral and deconfinement transitions [TJ [21 IT0rfT5] . Additional effects of the quark fields on the Polyakov loop potential 
and the phase transition at finite density have also been investigated [3H5] ■ 

A better understanding of the mechanism at the origin of these transitions requires the investigation of fluctuations 
in the PNJL model. In this paper we will show how fluctuations can be included by performing numerical simulations 
of the thermodynamics using standard Monte-Carlo (MC) techniques. The advantage of this method is that it 
automatically incorporates fluctuations to all orders. In the present work we restrict ourselves to fluctuations of the 
static zero-modes which lead to an improvement beyond the saddle-point approximation. 

Since we choose a finite volume for the Monte-Carlo evaluation, in this respect our approach is similar to lattice 
QCD calculations. Fluctuations of particular interest are those involving Goldstone modes of both zero and finite 
momentum. In principle, it is necessary to take all Goldstone fluctuations into account. These fluctuations will 
restore the chiral symmetry in a finite volume in the absence of explicit symmetry breaking. Such effects have been 
investigated for example by using Renormalization Group methods |16H19j . For sufficiently large explicit breaking of 
the chiral symmetry through a non-zero quark mass, as in the present work, it is legitimate to restrict the Goldstone 
fluctuations to the zero-modes. Finite-volume effects in the NJL-model have also been investigated in mean-field 
calculations [20-22 and using sophisticated Dyson-Schwingcr methods for QCD [2"5] . 

We perform our analysis in particular for the case of vanishing chemical potential where a comparison with lattice 
simulation results is possible. In this case we do not need to include diquark degrees of freedom in the model. The 
behavior of the chiral condensate and the Polyakov-loop expectation value are only slightly affected by the presence 
of fluctuations but their effect on the susceptibilities is much more pronounced. In particular, the temperature 
dependence of the flavor non-diagonal second derivative of the thermodynamic grand canonical partition function 
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with respect to quark chemical potentials is particularly sensitive to these fluctuations. For example, a quasi-particle 
model calculation [24 finds a vanishing result for this susceptibility. The same result is found in a PNJL calculation 
using the saddle point approximation [35]. An exception to this general behavior is discussed in the NJL approach 
of ref. |26| where a non-zero off-diagonal susceptibility is observed in the presence of an isovector vector coupling 
between quarks. In order to account for these flavor non-diagonal susceptibilities, we include a bosonic field with 
the quantum numbers of the pion. In the mean field approximation, the expectation value of this field corresponds 
to a pion condensate. Although this pion condensate vanishes at (i = 0, fluctuations of this field contribute to the 
flavor non-diagonal Taylor expansion coefficients for the susceptibilities. These zero-mode fluctuations appear only 
in the finite- volume system and vanish in the infinite- volume limit. In addition, fluctuations in the gauge fields 
also contribute to the non-diagonal susceptibilities. These fluctuations relate to the difference of the Polyakov loop 
expectation value and its conjugate at finite chemical potential [25] and persist even in the infinite volume limit. 
Gauge field fluctuations alone are not sufficient to explain current observations on the lattice. A combination of 
pionic and gauge field fluctuations, on the other hand, appears to succeed. 

This paper is organized as follows: the partition function of the PNJL model is introduced in Section|ll] In Scction [nT| 
we discuss the model in a finite volume. An additional parameter, the ratio between the spatial and temporal extent 
of the Euclidean volume, is introduced, analogous to the one in finite-temperature lattice simulations. In this context 



we apply the Monte-Carlo method to the PNJL model. In Section IV results for the pure gauge case are presented, 
using the Monte-Carlo method with a suitable Polyakov-loop effective potential. In Section [V] we study the behavior 
of the chiral condensate, the chiral susceptibility and the Polyakov loop for Nf = 2 quark flavors. Results concerning 



the Taylor expansion coefficients of the thermodynamic potential are shown and discussed in Section VI Conclusions 



are presented in Section VII 



II. THE PNJL PARTITION FUNCTION 

The Euclidean action of the two-flavor PNJL model including finite baryon and isospin chemical potentials is given 
by [HI HZ] 

S E (i>,M = f dT J d 3 x{^(^ + 7 ^-m)^ + G[(^) 2 + (^ 75 T^) 2 ]}-/3 J d 3 x U{<t>, 0), (1) 

with (3 — 1/T. Here ip is the Nf — 2 doublet quark field, m = diag(m u ,m<j) is the quark mass matrix and the 
covariant derivative is 

ip = i^(d»-igAi>). (2) 
The quark chemical potential matrix jl is defined as 



/' 



[id 



The Polyakov-loop effective potential U involves the gauge field degrees of freedom denoted by cf> and models the 
confincmcnt-dcconfincment transition in the pure gauge theory at mean-field level. In the PNJL model quarks 
interact with a background color gauge field A 4 = iA , where A = 8^gA^t a with the gluon fields A£ £ SU(3) C and 
t a = X a /2. The field A 4 is related to the traced Polyakov loop according to 

$ = -^-^ C L with L = exp (i f£ &tA±j . (3) 

In the Polyakov gauge, the matrix L is given in a diagonal representation 

L = cxp(i(0 3 A 3 + (Ms), ( 4 ) 

with the (diagonal) 5/7(3) generators A 3 and A 8 . The dimensionless effective fields 4> 3 and 4> s are identified with 
the Euclidean gauge fields in temporal direction divided by the temperature, A±' /T and A± /T. These two fields 
parametrize the diagonal elements of SU(3) C . 

Such a description of the pure gauge thermodynamics is supposed to be valid in a limited range of T. It is known 
that, at very high temperatures, the explicit presence of transverse gluons begins to govern the dynamics. As a 
consequence, the validity of our approach is limited to temperatures less than about 2 T c , where T c ~ 0.2 GeV is the 
typical transition temperature scale. 
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In this paper we consider the ansatz for the effective potential given in motivated by the SU(3) Haar measure 

which happens when integrating out six of the eight gluon fields: 

U($,§*,T) = _ l a ( T ) $ , $ + ^ ^ _ g $ * $ + 4 ( $ *3 + $ 3) _ 3 ($*$)2] ( 5 ) 

The temperature-dependent prefactors are given by 

a(T) = a + ai + a 2 and b(T) = b 3 . (6) 



As will be shown in Section IV the particular choice of a(T) and b(T) is such that we can reproduce the high- 
temperature behavior of thermodynamic quantities like pressure, energy and entropy density. In this way we indirectly 
take into account effects connected with gluonic degrees of freedom that are integrated out in the definition of our 
effective potential. An additional constraint for fixing the parameters is the critical temperature of the first-order 
deconfinement transition in pure gauge QCD, T = 270 MeV, as given by lattice calculations, and the requirement 
that $*,$-> 1 as T-> oo. 

Given the action ([I]), the partition function of our system is 

Z = M fvmV^exp(-S E [^,iJ,^) , (7) 



where <j) stands for the Polyakov loop fields 03 and cf> & . Applying standard bosonization techniques, multiplying the 
partition function by the expression 



T>aT>TT exp 



d 



ml m 



and evaluating the resulting Gaussian integral over the fermionic degrees of freedom, one finds 

Z = M J V$DoVk dettS" 1 ] exp [ - ^ / d 3 x(u(<b, T) + ^ ^ )] ■ (9) 
We write the pion field tt = (ir 1 , tt 2 , tt 3 ) in terms of tt^ 1 = -^{n 1 ± Z7r 2 ), tt° = tt 3 and r ± = ^(t 1 ± ir 2 ), so that 

f • 7T = V2(t+7T~ + T~7T+) + T 3 7T°. 

The inverse quark propagator takes the form 

g-i = ( ~$ + (A*™ - + i757r° - M iV2<y 5 %+ \ 

\ iV§Tto*- -Q + - iA 4 )j Q + i 75 7r — M J 

with the dynamical quark mass M = ttiq — a generated by the scalar field a < 0. We work in the isospin symmetric 
limit with m — m u = for convenience. This scalar field is related to the chiral (quark) condensate by a — G(ipip). 
Equivalently, the partition function can be rewritten as 

Z = N J P0P<7X>7rcxp(-S[cr,7r,A 4 ]) 

= N J V^VaVnexp ^Trln^" 1 ] - ± / d 3 x(u(^T) + ^ !L ) • (10) 

The standard path for the evaluation of this partition function is to first consider the mean-field limit, where only 
the classical trajectory in the path integral is taken into account. In this approximation, the partition function 
becomes an ordinary multi-dimensional integral. This mean-field partition function is evaluated using the saddle- 
point approximation, taking into account only the maximum contribution to the partition function. To evaluate this 
contribution, the action is minimized by variation of the classical fields. This approximation becomes exact only in 
the infinite-volume limit, which is not reached in lattice simulations. In the present paper we will consider a different 
approach. 
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III. PNJL MODEL IN A FINITE VOLUME 
A. Finite- Volume partition function 

In the present calculation we perform a step beyond mean-field approximation by including fluctuations of the 
zero modes of the relevant fields. This is admittedly only part of all possible field fluctuations, but it represents 
nevertheless an improvement with respect to the standard mean-field calculation. These zero-mode fluctuations can 
be introduced considering a system defined in a finite volume V. 

The partition function in momentum space is written as 

Z= /^^expf^^^Trln^^^-W^T)-^ 2 ^)] (11) 

n p 

where u> n = (2n + l)nT are the Matsubara frequencies. Given the inverse quark propagator defined by 

g-lfj^ $\ = ( ( iuj n + M« - ^4)70 + ij5^° - 7 -P- M iV2j 5 TT+ 

V iV2j 5 n' (iiu n + n d - iA 4 )j + ij 5 TT° - j ■ p - M 

the fermionic determinant can be evaluated by diagonahzation and one finds 

ilHn[S-] = -2AT, J |^£{Tln[l + e-^] + \aE 3 } (12) 



with twelve quasiparticle energies 

£1,2 = \fW) + M/) 2 + ± (ji - 2iA 8 /V3), E 3A = yj{e{p) - in) 2 + 2n+n- ±(p- 2iA s /V3), 

E 5S = y/( £ (p) + tn) 2 + 2tt+tt- ± (fi + i{A 3 + A s /V3)), E 7 , s = y/(e(p) - in) 2 + 2n+n- ±(p + i{A 3 + A 8 /V3j), 

£9,10 = V(£(P) + A*/) 2 + 27T+7T- ± (/x - i(A 3 - As/VS)), £11,12 = y/(e(p) - /i/) 2 + 2it+tt- ±(ji- i{A 3 - A 8 /V3j), 

where we have introduced isoscalar and isovector chemical potentials, \i = (/i n + /irf)/2, /i/ = (p u — fJ.d)/2 and 
e(p) = \/ p 2 + 7Tq + (mo — a) 2 . The energy difference is defined as the difference between the quasiparticle energy and 
the energy of a free fermion, £0 = \/p 2 + ± /i/: AEj = Ej — Eq =F fj,. 

For the volumes considered in the present work, effects from discretization of momenta turn out to be negligible. 
In order to simplify the numerical evaluation of the partition function, we have therefore substituted an integral for 
the sum over the discrete three-momenta. This approximation becomes invalid in the limit of small volume size and 
small quark masses and it is not applicable in general. 



The presence of a volume factor V in the exponent of Eq. (Ill makes it possible to compute the full partition 
function in mean-field approximation using Monte-Carlo techniques, as explained in the following section. In this way 
we consider not only the saddle-point contributions, but also configurations that correspond to fluctuations around 
the minima of the action. The zero-mode fluctuations of the pion fields which appear at this level will disappear in the 
infinite-volume limit. Additionally, contributions from the gauge fields which are related to the difference between the 
Polyakov loop and its conjugate also contribute to flavor non-diagonal susceptibilities. These contributions survive in 
the infinite- volume limit. 

The size of the volume is now specified according to the conventions adopted in lattice calculations. For a fixed 
extension of the lattice in the Euclidean time direction, the temperature is set by the lattice spacing a, and the volume 
size is related to the temperature: 

1 N 3 

a =wf v = N ° a3 = w^> (13) 

where N t is the number of lattice sites in the Euclidean time direction, and N s is the number of lattice sites in the 
space direction. It follows that 

V = k/T 3 , (14) 

where different values of k = (N s /N t ) 3 will be chosen for our purpose. On the lattice the parameter k must be 
sufficiently large in order to reduce finite volume effects while keeping the system at a specified temperature. On 
the other hand, for increasing values of k the necessary computing time grows exponentially at fixed temperature 
and simulations for very large volume become less feasible. At present, N s /N t = 4 (k = 64) is the largest value 
used in lattice simulations [29] [30] . In the MC-simulation of the PNJL model, the aim is to understand the volume 
dependence of thermodynamic quantities. For this reason we perform calculations at different values of fc, including 
k = 64 which corresponds to the typical lattice simulation volume. 
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B. Monte-Carlo method 



In Euclidean quantum field theory, the expectation value of an observable O is given by 

J T>(pe- S M 

where ip stands for the set of relevant field variables. A statistical method is used in order to select those configurations 



of the fields tp which give a significant contribution to ( 15 ). An efficient way of computing the ensemble average consists 
of generating a sequence of configurations with a probability distribution given by exp(— <S[y]). This technique is called 
"importance sampling" . If the sequence generated constitutes a representative set of configurations, then the ensemble 
average (O) will be given approximately by the sum 

1 N 

(°)^^° {(pi) ' (16) 
»=i 

where (pi (i = 1, ...,N) denote the configurations generated. 

There are different techniques to generate the relevant ensemble of configurations. In this paper we adopt the 
Metropolis method. The algorithm is simple: starting from a random configuration C, a new configuration C' is 
generated with acceptance 

p = min{l,exp(-5[C"])/exp(-5[C])}. (17) 

The new configuration will be accepted if the action has decreased; if it has increased, it will be accepted with 
probability exp(— S\C'})/ exp(— <S[C]). If we use a cooling procedure and accept only those configurations that decrease 
the action, the chain of configurations generated in this way ends up with a configuration that corresponds to the 
saddle point result, i.e. with a configuration that minimizes the action. In the next sections we will see how this 
Monte-Carlo evaluation of the PNJL model works both in the pure gauge sector - i.e. with the Polyakov loop potential 
- and in the case with Nf — 2 quark flavors. 



IV. MONTE-CARLO EVALUATION OF THE POLYAKOV LOOP SECTOR 



The effective potential for the Polyakov loop in Eq. ((si) is written as a function of the traced Polyakov loop $. In 

A <3} A i8} 

this section we study the properties of this potential in terms of the fundamental gauge fields 4>3 = —jr- and 0s = . 
Using the definition of the Polyakov loop given in Eq. ([3]), and the representation Q 



we obtain for the effective potential ([5]) expressed in terms of 3 and <fig: 

^ (03 r t 8 ' T) = -^a(T)[3 + 2cos(20 3 )+4cos(^)cos(V30 8 ) 



+b(T) In 



16 
27 



(cos(0 3 ) - cos(V3</> 8 )) sm((f> 3 )' 



(18) 



(19) 



In this section the Monte-Carlo approach is applied to the sampling of the Polyakov loop. Although contributions 
from fluctuations are small in this case, it is instructive to illustrate how the approach works. We start from a generic 
configuration fixing the A\ and A 4 fields at a given temperature. Using the Monte-Carlo method, we generate 
points in configuration space close to the minimum of the action. In particular, applying cooling, we reach exactly 
the minimum. Fig. [I] illustrates the difference between the trajectories generated by the cooling procedure and the 
field samples produced by the Metropolis algorithm at the critical temperature T = T = 0.27 GeV. For cooling (left 
figure) the action must decrease, and in a few steps one of the three minima of the potential is reached. On the other 
hand, using the full Metropolis algorithm, the action can also increase occasionally, and as a consequence the space 
of the available configurations is much larger (right figure). 

To fix the parameters of the effective potential, we have computed the pressure, energy and entropy density for the 
pure gauge theory at different temperatures, including fluctuations, and chosen parameter values to fit the available 
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FIG. 1: Comparison between the trajectory generated by a cooling algorithm (left) and the field configurations produced by 
the Metropolis algorithm (right) applied to the Polyakov loop effective potential at temperature T = T c — 0.27 GeV. The lines 
represent equipotential curves of U (Eq. |5|) in the ^-^Is-plane. 

lattice data (Fig. |2] left). It is found that only the value of the parameter 63 differs slightly from the one obtained by 
mean-field calculations in the saddle point approximation. With the smallest value for the parameter k used in the 
calculation, which leads to the largest fluctuations, we find 63 = —1.64, compared to 63 = —1.68 from the saddle point 
approximation. In the right panel of Fig. [2] we compare the expectation value for the Polyakov loop obtained in the 
saddle point approximation (line) with our results (open points) and the lattice data (points). There is practically 
no difference between the two methods, cooling versus Metropolis: both agree well with lattice data. Once we have 




FIG. 2: Left: Results for pressure p, energy density e and entropy density s from the Monte-Carlo evaluation of the Polyakov loop 
effective potential fitted to lattice results (lines), compared to lattice QCD data from [31] (open symbols). Right: Comparison 
of the Polyakov loop expectation value from the saddle point approximation (solid line) and the Monte-Carlo evaluation (open 
points) with lattice QCD data taken from [32] (solid points). 

established how the Monte-Carlo approach works in the case of Polyakov loop dynamics, we can move on to the more 
interesting case with two quark flavors. 



7 



V. MONTE-CARLO APPROACH TO THE TWO-FLAVOR PNJL MODEL 



The starting point for studying the thermodynamics for Nf — 2 quark flavors is the partition function (111. The 



degrees of freedom in this case are the A^ 3 -* and components of the gauge field, and the bosonic field variables a 
and 7? . In th e ev aluation of the path integral, we need to fix the volume V as a function of the temperature. Looking 



back at Eq. (14 1, we see that this requires fixing the dimensionless index LT = fc 1 / 3 of the Euclidean volume. In this 
work we consider six different choices of fc. The value fc = 64 corresponds to the largest one currently used in lattice 
simulations. In addition we consider fc = 125, 250, 500, 1000, 2500. The ratio between the smallest and the largest 
fc is ~ 40. In this way we can study systematically the dependence of the observables on the volume size at fixed 
temperature T. In the NJL sector of the model we also need to specify the current quark mass mo, the coupling 
constant G and the three-momentum cut-off A. The parameters used here are the ones of Refs. [T2l IT4"] : 

m = 5.5 MeV, G = 10.1 GeV~\ A = 650 MeV. 



A. Chiral and deconfinement transitions 



The principal aim of this work is to contribute to the investigation of the QCD phase diagram. In this section we 
study how the chiral and deconfinement transitions arc affected by the introduction of fluctuations around the mean 
field in our Monte-Carlo PNJL approach. This is achieved by evaluating the chiral condensate and the Polyakov 
loop expectation value for different volumes and comparing with the mean-field result in saddle point approximation. 
This comparison is presented in Fig. [3] The presence of fluctuations does obviously not modify the behavior of the 
Polyakov loop expectation value; the four different sets of data overlap perfectly. For the chiral condensate below 
the critical temperature, we notice that there is in fact a non-trivial dependence on the temperature: the expectation 
value of the field a starts to decreases earlier for smaller fc. This reflects a volume dependence that moves the chiral 
transition temperature from T a — 254 MeV for fc = 2500 to T a = 222 MeV for fc = 64, as deduced from the following 
analysis of the chiral susceptibility. 




FIG. 3: Dependence of the chiral condensate, (t/>-!/j)t/(t/>7/>}o, and of the Polyakov loop expectation value (<&) on the parameter 
k — (LT) 3 in a finite volume. Deviations from the mean-field result in the infinite-volume limit are manifest only for the chiral 
condensate. The behavior of the Polyakov loop is completely unchanged. 



B. Chiral susceptibility 



The chiral susceptibility is sensitive to fluctuations in the fields and therefore provides a good test for the Monte- 
Carlo evaluation of the PNJL model. In the infinite- volume limit, V — > oo, the Monte-Carlo calculation should recover 
the saddle-point result. To perform this comparison, we also calculate the chiral susceptibility in the saddle-point 
approximation. 



The chiral susceptibility is denned as 



Xa 



T d 2 
V dniQ 



In Z(m ,T), 



(20) 



in terms of the partition function Z of Eqns. ( 10 ),( 11 1. The second derivative is taken with respect to the quark mass 
m . Performing these derivatives leads to 



Xa = 



V 
T 



+ 

V 
T 



1 

Z(m ,T) 
1 

Z(m ,T) 
1 

Z(m ,T) 



VaVnVA 



VaVnVA 



dlndet S r - 1 (m Q ,T, a, n,A) 
dniQ 

d In dct S~ 1 {m Q , T, a, if , A) - S [^,A] 



-S[<t,k,A] 



dmo 



VaVirVA 



d 2 lndetS , - 1 (m ,7 1 , a, w, A) 



-5[ct,tF,A] 



dlndet S~ x {mQ,T, a,T?, A) 
dm 

d 2 lndet S , - 1 (m ,T, a, tt, A) 
dm 2 , 



d In det S~ 1 (mo ,T,ej,n,A)\ 2 ' 
dm I . 



(21) 



This expression can now be evaluated using the Monte-Carlo algorithm. We use the position of the peak in the chiral 
susceptibility as a measure for the chiral transition temperature T a . The results presented in Fig. [4] show that T CT 
moves toward its saddle point limit, T a w 254 MeV (black points). 

The mean-field numerical result is obtained by evaluating the thermodynamic potential in saddle-point approxima- 
tion for current quark masses mo in the range from 1 to 10 MeV and for different temperatures. The thermodynamic 
potential is then interpolated and the second derivative is calculated numerically. The result is represented by the 
solid curve in Fig. [5] The Monte-Carlo calculation approaches the mean-field limit when the volume is increased at 
fixed temperature. 
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FIG. 4: The chiral susceptibility as a function of temperature for different volume aspect ratios. We use the peak position as 
a measure for the chiral transition temperature T„. The solid line shows the trajectory of the chiral transition temperature as 
it rises with increasing volume. 



The role of fluctuations in the Monte-Carlo calculation of the chiral susceptibility can be understood from Eq. (21 1. 
The disconnected contribution (the term in the square brackets), vanishes in the infinite- volume limit since fluctuations 
of the mean field contribute as ~ 1/V; 

(• 2 ) - (») 2 ^0 as V -> oo, (22) 



Because the prefactor V/T in Eq. (21 1 compensates the volume dependence of the leading fluctuation contributions, a 
finite susceptibility results also in the limit V —¥ oo. Since additional contributions of fluctuations in the mean fields 
are of higher order in 1/V, the saddle point approximation becomes exact in this limit. 
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X, 




° 0.8 



T/Ta 

FIG. 5: Chiral susceptibility as a function of T/T a : With increasing volume the Monte-Carlo results (open symbols) approach 
the saddle-point mean field result (solid curve). 

VI. NON-ZERO QUARK CHEMICAL POTENTIALS: TAYLOR EXPANSION 

Dealing with non-zero quark chemical potentials /j, q in lattice QCD thermodynamics is notoriously difficult because 
of the well-known fermion sign problem. A possible way of overcoming this problem is the Taylor-expansion approach. 
Instead of performing an explicit calculation at [i q ^ 0, the thermodynamic potential is expanded in a Taylor series 
in H q /T around zero chemical potential, 



n(r,/i) 



i 



VT 3 



lnZ=£ 



with 



Xij{T) = 



i,j=0 



d i+j n 



T ) 



(23) 



(24) 



where only even terms survive due to CP symmetry. The coefficients XijiT) are evaluated at [i q = 0. 

The comparison between lattice data and Monte-Carlo calculations for these coefficients in the PNJL model represent 
an important test of this model. In particular the flavor non-diagonal coefficient \n that vanishes in the saddle point 
approximation is of interest in this context: it is necessary to take fluctuations of the mean field into account in order 
to obtain a non- vanishing result for xn- Since the strength of fluctuations depends on the volume, we again evaluate 
the Taylor coefficients for different volume sizes at each temperature, i.e. for different values of the parameter k. 



A. Second order Taylor expansion coefficients and susceptibilities 



The first derivative of the logarithm of the partition function gives 



dlnZ(T,/i u ,fi d ) 



' ' In f VaVirVA exp lndct 5 _1 (T, fj, u , /j, d , cr, if, A) — S g [A] 



Z(T,fi. 

V / d lndct g'^T, fi u ,fi d ,(T,Tv,A) \ 
T\ dn q " ~ I 



A] 



(25) 
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Proceeding in the same way for the second derivative, we obtain the coefficients Xuu and Xud (quark susceptibilities) 

*" = ^a^ lnZ = ^(H5^ lndet ^ 1(T ' Mu ' Md '"' 7f ' j4) 



T J \\ dfi u J \T J \ dfj, u 

Consider first the flavor-diagonal susceptibility X20 — \Xuu- A comparison between Monte-Carlo results and the 
saddle point approximation for Xuu is shown in Fig. [6] for different volumes as a function of the temperature. In this 
case the contributions of fluctuations arc evidently small, reflecting the fact that Xuu is governed by the non-vanishing 
quark condensate and the A 3 component of the gauge field. 

The behaviour of the flavor non-diagonal coefficient xu — Xud, on the other hand, is quite different. It vanishes in 
the saddle point approximation whereas lattice QCD clearly displays a non-zero signal for this quantity around T c . 




FIG. 6: Flavor-diagonal second order expansion coefficient x«« = 2^20 for different ratios k compared with the saddle point 
approximation. The Stefan-Boltzmann limit is marked by SB. 



For the detailed evaluation of Xud we return to the expressions, Eqs. (25), (26 1, and evaluate derivatives of the 
fermion determinant: 



and 



d 
dp u 
. 1 



bidetS 1 {T, fin, fid, <r,n,A) 
sin Ai 



dp p 



sin Ao 



sin A^ 



cos Ai + cosh(e (p) / T) cos A 2 + cosh(e (p) / T) cos A 3 + cosh(e (p) / T) 



d 2 
dpudpd 



IndetS 1 (T, p u ,p dl a, if, A) 



1 + cos Xi (cosh(e(p) /T) cos Aj 

Te 2 (p)(cos A, + cosh(e(p)/T)) 2 + e 3 (p)(cos A, + cosh(e(p)/T)) 



cos 

"73 



h( £ (p)/T)-sinh( £ (p)/r) Y 



dp 



3p 2 



»(p)(coBA i + cosh(e(p)/T))yy tt 2 J f e{pf 
where e(p) = \J p 1 + (mo — <r) 2 + 7? 2 , while the A^ are the eigenvalues of the Polyakov loop matrix: 



. _ ^3 . A s A 3 A 8 

A l — ~ H 7^,, *2 — H ^3 



T y/3T' 



T J3T 



2A S 



(27) 



(28) 



(29) 
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From Eq. ( 27 1 we observe that this expression is odd with respect to As and even in all other fields. As a consequence 
the expectation value of such a term is zero when the functional integration on the domain of Ag is performed. This 
implies that the third term in Eq. ( 26 ) vanishes. 



Another crucial observation is that taking the mean field limit for the pion field, it = (tt) = 0, Eq. (28) vanishes. 
Consequently, Xud vanishes altogether in the mean field approximation, independent of the temperature. Computing 
the expectation values of Eqs. (27) and (p8|) using the MC-PNJL approach we include corrections induced by fluctu- 



ations of the pionic and Polyakov loop fields. Moreover, the main contribution to the Eq. ( 28 1 is given by the pionic 
fluctuations, whereas the Eq. (27) is non-zero mostly due to fluctuations of A$. 

The pionic and As contributions to Xud resulting from the MC-PNJL computation are shown in Fig. [7| Two 
characteristic features are immediately apparent. First, the term involving pionic zero-modes is strongly volume 
dependent and vanishes in the limit of infinite volume. Secondly, the term associated with fluctuations of the A% 
gauge field is independent of the box size and survives in fact as the volume becomes infinitely large. 




FIG. 7: Different contributions to the off-diagonal susceptibility xn — X«<* f° r different volume ratios k computed in the 
Monte-Carlo approach. Left panel: Contribution from pionic fluctuations, for which the volume dependence is large. Right 
panel: Contribution from fluctuations of the As field, which show a negligible volume dependence. 



B. Chiral effective Lagrangian 



In order to better understand the role of the pionic fluctuations in the evaluation of Xud, let us briefly digress and 
study this issue in the context of chiral perturbation theory (ChPT). 

For low temperatures and small values of the chemical potential, the physics is dominated by the effects of light 
pions and we can describe the system in terms of an effective chiral Lagrangian. While not influenced by the baryon 
number chemical potential, the pions do couple to the isovector chemical potential and its effect can be included in 
this effective Lagrangian. The chemical potential enters into the QCD Lagrangian like the zeroth component of a 
gauge field [531 133- When one promotes the global chiral flavor symmetry SU(2) L x SU{2) R of the QCD Lagrangian 
to a local symmetry, gauge invariance determines completely how the chemical potential must be implemented in the 
effective chiral Lagrangian [35H37j . This Lagrangian has the form, expressed in terms of the chiral field £, 

C = ^Tr [V.EV,^] - ^Tr [£ + £+] , (30) 

where £ = | (ipip) | is the magnitude of the chiral condensate and the covariant derivatives are defined as 

V E = 9 E + m [r 3 E - Er 3 ] , V,E = d. t E, i = l,2,3 

v e+ = a Et + m r 3 Et - Etr 3 ] , v 4 E f = «%£+, i = 1,2, 3. (31) 



T3 = diag(l, —1) is the diagonal isospin generator. Inserting these expressions into the chiral effective Lagrangian (30) 
allows to identify the terms depending on the isovector chemical potential: 

C = ^Tr [c^E^Et] +^ M/ Tr[(a E)EtT3 + Et(a E)T 3 ] 

+ | M7 Tr (r 3 ET 3 Et - 1 2 ) - ^Tr [E + £t] (32 ) 
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For > m w /2, the formation of a pion condensate is possible and the ground state changes from the one at /i/ = 0, 
which is determined solely by the chiral condensate, to a combination of chiral and pion condensates [35, 36j 138) . 
Since we are interested in the determination of the susceptibility \ud defined in the limit fii = 0, we can expand 
around the unchanged ground state and the chiral field E £ SU(2) can be parameterized as 



(33) 



E = exp [ &-7-T* 

Jit 



where the r a , a = 1, 2, 3 are the generators of SU(2) with normalization Tr(r a r b ) = 25 ab . Expanding the Lagrangian 
to second order in the pion fields 7r a , one finds 



£ 



1 



2„2\ 



--m 2 7r a 7T a - 2/4(ttV + tHtT 



(34) 



where we have identified toE = f 2 m 2 . In order to make contact with the results from the Monte-Carlo evaluation of 
the PNJL model with fluctuations of the mean-field only, it is sufficient to take the static part of the Lagrangian into 
account. The partition function for this static part is 



3 

f f[dTT a exp< 


r v 


' a=l 


3 

f X\d-K a exp< 

' a=l 


r v 



m^Tr" - 2/i?(7T 1 7T 1 + TtV) 



The second-order off-diagonal expansion coefficient %ud is given by 

1 d 2 



(tt) 



lnZ s 



(35) 



(36) 



and receives contributions from the two-point correlations (tt + tt ) . Evaluating the integral, we obtain the static result 
f° r Xud from the chiral effective Lagrangian: 



2 

VTml 



2 T 2 1 

k m 2 



(37) 



setting again V — k/T 3 . This prediction can be compared directly with our Monte-Carlo PNJL results, provided we 
take the temperature dependence of the pion mass into account, using the relation given in |39j : 



J7v(T) =m n [ 1+ — 2 YO{y 



9i(ml,T,L)-- ^ 
n = (niL, n2L, n^L) 



— / dA'- 3 ^exp(-m 2 A 



n 2 /(4A)), 



(38) 



(39) 



Fig. [8] demonstrates that the picture so obtained from the chiral effective Lagrangian is completely consistent with 
our Monte-Carlo calculations in the PNJL model, as far as the pionic contributions to \ud are concerned. From Fig. [8] 
it also follows that the chiral perturbation theory prediction for this coefficient is reliable until around T/T c ~ 0.7. 



C. Comparison with lattice data 

Lattice QCD studies of Xud have been carried out for example in Refs. [25II3D], both with k = 64 but with different 
quark masses, corresponding to pion masses m, = 230 MeV and 770 MeV. These lattice results are compared to 
our Monte-Carlo PNJL computations (using the physical pion mass) in Fig. [9] The shape of the \ud signal is quite 
well reproduced within the large error band of the lattice data. The difference between lattice results computed with 
different pion masses is now quite plausible. Given that the pionic fluctuations dominate over those from the A$ 
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component of the Polyakov loop field, this behavior is just what one expects from Eq. (37). At the same time one 



would expect that lattice simulations performed ideally with physical quark masses would actually yield even larger 
magnitudes of Xud than those with m w = 230 MeV. The Monte-Carlo results notably include only the pionic zero 
modes. Finite-momentum fluctuations would tend to further increase the pionic effects in Xud- 

In the infinite- volume limit, only the gauge field signal (plus possible finite-momentum pionic modes) survives in 
Xud, as pointed out earlier. It would be interesting to see whether this predicted behavior is realized in lattice QCD 
when moving towards larger values of the ratio N s /N t . 



MC-PNJL 
m n =0.139 GeV 




& Lattice 

m n =0.770 GeV 

■ Lattice 

mn=0.230 Gev 



T/Tc 



FIG. 9: Temperature dependence of the flavor off-diagonal susceptibility Xud in the Monte-Carlo approach to the PNJL model, 
using k — 64 (LT = 4). Lattice data [291 130| with the same volume aspect ratio LT and different pion masses are also shown 
for orientation. 



VII. CONCLUSIONS 



Spontaneous chiral symmetry breaking and confinement are well-known phenomena emerging in the study of QCD 
thermodynamics. Both these features are incorporated in the Polyakov-loop extended Nambu-Jona-Lasinio (PNJL) 
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model. Predictions for the deconfinement and the chiral restoration transition temperatures, obtained from mean- 
field calculations in this framework, compare quite well with the available lattice data. Nevertheless, the mean-field 
approximation is not sufficient if one wants to correctly account for the temperature dependence of other observables 
such as the Taylor expansion coefficients of the pressure in powers of u-and d-quark chemical potentials. 

In this work we have applied standard Monte-Carlo techniques to a PNJL model in order to go beyond the saddle- 
point approximation. This becomes important when the system is considered in a finite volume. The strength of the 
fluctuations introduced in this way depends on the size of the volume. At fixed temperature T, the spatial volume size 
L = V 1 / 3 is determined by the factor fc 1//3 = LT of the Euclidean volume. We have checked the method by studying 
the chiral susceptibility for different volume sizes V and shown that the saddle-point result is approached in the limit 
V — > oo. Studying the thermodynamics, we found that the introduction of mean- field fluctuations in a finite volume 
does not change the traced Polyakov loop, both in the pure gauge sector and for Nf = 2 quark flavors. The chiral 
condensate with Nf = 2 experiences modest changes through such effects. 

On the other hand, the inclusion of such beyond mean field fluctuations in a finite volume does affect the susceptibil- 
ities significantly. We find that their impact is crucial for the evaluation of higher-order Taylor expansion coefficients 
of the pressure. In particular, the second-order flavor non-diagonal expansion coefficient Xud becomes non-zero. Our 
result from a Monte-Carlo computation agrees well with lattice data using the same k for the Euclidean volume. In 
contrast, this coefficient vanishes in the saddle-point approximation. Good agreement with lattice results is also found 
for the second diagonal moment of the pressure. These results show that finite- volume fluctuations with pion quantum 
numbers can account for the non- vanishing expectation value of the susceptibility Xud- This zero-mode contribution 
vanishes in the infinite- volume limit. Concerning the role of the pionic zero- mode fluctuations, it is also demonstrated 
that the Monte-Carlo results for Xud are fully consistent with those from chiral perturbation theory for temperatures 
below T c . 

The contribution from the difference in the expectation values of the Polyakov loop and its conjugate, which 
also contributes to Xud through the fluctuations of the Ag component of the gauge field, remains finite even in the 
limit V — > oo. In lattice simulations, both types of mean-field fluctuations (pionic and gauge fields) contribute to the 
observed susceptibilities, although the pionic effects are expected to be smaller for the larger pion masses used in lattice 
QCD. While these pionic fluctuations will be suppressed for large volumes (i.e. a large LT at fixed temperature), 
gauge field effects can still account for a non-zero Xud even in the infinite- volume limit. The present analysis predicts 
a decrease of Xud for larger volume sizes in lattice simulations. 

In conclusion, we have shown how finite-volume fluctuations in the PNJL model at non-zero chemical potential can 
be studied using a Monte-Carlo evaluation of the partition function, and that such fluctuations can affect suscepti- 
bilities significantly. Going beyond the mean-field approximation for the investigation of critical behavior requires 
additionally to take fluctuations on all momentum scales into account, and such a calculation would also be amenable 
to a Monte-Carlo evaluation. 
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